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ABSTRACT 

We present the structure of the 3D ideal MHD pulsar magnetosphere to a radius ten 
times that of the light cylinder, a distance about an order of magnitude larger than 
any previous such numerical treatment. Its overall structure exhibits a stable, smooth, 
well-defined undulating current sheet which approaches the kinematic split monopole 
solution of Bogovalov 1999 only after a careful introduction of diffusivity even in 
the highest resolution simulations. It also exhibits an intriguing spiral region at the 
crossing of two zero charge surfaces on the current sheet, which shows a destabilizing 
behavior more prominent in higher resolution simulations. We discuss the possibility 
that this region is physically (and not numerically) unstable. Finally, we present the 
spiral pulsar antenna radiation pattern. 

Key words: pulsars: general — Magnetic fields — radiation mechanisms: general — 
Methods: numerical — Magnetohydrodynamics (MHD) 



1 INTRODUCTION 

Almost half a century since the discovery of pulsars and de- 
spite the great observational progress, a comprehensive un- 
derstanding of the observed radiation remains elusive. This 
is to a large extent due to its breadth in frequency (from ra- 
dio to 7— rays) and dependence on the pulsar rotation phase 
as well as on its potential dependence on the (unknown) ob- 
server inclination angle and the angle between the rotation 
and magnetic axes. 

Pulsar radio emission remains the hardest to model 
(mainly because of its apparently coherent character) and 
most models that claim to account for its main features 
(multiple pulse profiles, polarization angle sweeps, pulse 
width-to-wavelength relation, etc.) are simply phenomeno- 
logical. It is considered that radio emission originates in 
(probably hollow) conical beams above the magnetic po- 
lar caps. The situation is slightly better at higher ener- 
gies (optical to gamma-rays) where the recent Fermi ob- 
servations indicate the emission region to be associated 
with t he outer magnetosphere even beyond the light cylin- 
der (R ansom et"alTl 20 id : [Guillemot et al.ll201lh . The exact 
emission mechanism remains unclear, however, it is generally 
considered to be curvature radiation emitted by relativistic 
electrons moving along the pulsar magnetic field lines near 
the separatrix between the closed and open field lines, al- 
though synchrotron and Inverse Compton emission could 
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potentially make important contributions to the spectrum 
(jMuslimov k Harding II2004I : iHarding et al.ll2008l ). 

Given the relativistic unidirectional motion of the pul- 
sar radiating particles, any significant advance in under- 
standing their emission across the electromagnetic spectrum 
must take into account the detailed structure of the pul- 
sar magnetosphere. Its structure, even in the axisymmetric, 
aligned dipole field had been unknown until rather recently. 
Following the original solution of Contopoulos, Kazanas & 
Fendt (1999), a number of additional treatments have es- 
tablished the precise structure of this configuration (Gruzi- 
nov 2005; Contopoulos 2005; Timokhin 2006). More recently 
the structure of the non-aligned dipole was obtained nu- 
merically (Spitkovsky 2006; Kalapotharakos & Contopou- 
los 2009); this was then used to produce kinematic model 
light curves by placing the emission of radiation both near to 
(Bai & Spitkovsky 2010) and well outside the light cylinder 
(Cont opoulos Ka lapotharakos 1 12010[ ). with results that 
apparently are both consistent with observations. The idea 
that the observed radiation is sited outside the light cylinder 
was not new. Lyubarsky (1996) had proposed the possibility 
that the high energy emission comes from sites just beyond 
the light cylinder at the reconnection region near the Y- 
point while Petri & Kirk (2005) proposed the stripped wind 
model in which the observed radiation comes from a zone 
far beyond the light cylinder. 

The most recent multi-frequency observations of pulsar 
emission do not seem to provide a preference between these 
rather diverse sites of pulsar emission, being in general phe- 
nomenological agreement with both. In fact, in most pulsars 
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with both gamma-ray and radio observations available (56 
as of April 2011), the radio pulse precedes the m ain ^amma- 
ray p ulsar by about 15% to 40% of the period (lAbdo et alJ 
I2OIOI ). This is consistent with the radio emission being as- 
sociated with the polar cap, and the high energy emission 
being associated with a region further out along the open 
field lines. Still, 7 pulsars (the Crab and another 6 out of 
20 known millisecond gamma-ray pulsars) do not follow this 
simple picture. In these, the radio and high energy peaks are 
in phase. This suggests that, at least in these 7 pulsars, ra- 
dio emission is associated with the extended outer magneto- 
sphere. Interestingly enough, in the case of the Crab pulsar, 
a radio component associated with the polar cap emission 
region has also been observed, albeit only at low frequencies 
(Moffett & Hankins 1996). This may further suggest that, 
at least under certain circumstances, radio emission may be 
coming both from the polar cap and the outer magneto- 
sphere. 

This issue as well as the general stability of the mag- 
netosphere at larger distances motivate a study of its struc- 
ture at larger radii. Furthermore, a theoretical model that 
purports to account for the origin of radio emission in the 
outer magnetosphere was first proposed by Ardavan (1998) 
(see also Ardavan et al. 2008) . Ardavan proposed that a su- 
perluminally corotating pattern of electric charges and cur- 
rents generates at large distances regions of enhanced elec- 
tromagnetic radiation in vacuum where the local Poynting 
flux drops as the inverse [and not the inverse square) of the 
distance from the source. We were particularly intrigued by 
this proposition and decided to investigate its implications 
in the pulsar magnetosphere through extensive numerical 
simulations that probe its structure well outside the light 
cylinder. As we will see, our results, although still not 100% 
conclusive, do not seem to support the Ardavan model. 

The study of the outer magnetosphere (up to lOr^c; 
ric is the light cylinder radius) requires, on one hand high 
enough resolution simulations in order to guarantee minimal 
numerical dissipation, and on the other the introduction of 
some (small) amount of diffusivity in order to reduce noise 
and stabilize the numerical scheme. These simulations ex- 
hibit a stable and well-defined undulating smooth current 
sheet tending towards the split monopole solution (Bogov- 
alov 1999). However, the non-diffusive simulations exhibit a 
destabilizing spiral region at the crossing of two zero charge 
(null) surfaces that happen also to be located on the cur- 
rent sheet. This destabilizing behavior is enhanced in sim- 
ulations of increased resolution suggesting a physical (and 
not numerical) instability. 

This paper is structured as follows: in § 2 we discuss 
our numerical procedure in comparison with previous treat- 
ments of the problem. In § 3 we present the pulsar antenna 
radiation pattern. Finally, our conclusions are summarized 
in §4. 



2 NUMERICAL SIMULATIONS 
2.1 Procedure 

The 3D force-free structure of the pulsar magnetosphere in 
ideal MHD has been investigated numerically only very re- 
cently by two groups (Contopoulos & Kalapotharakos 2010, 



hereafter CKIO, and Bai & Spitkovsky 2010). Both inves- 
tigations have been based on the formalism of Force-Free 
Electrodynamics (hereafter FFE; Gruzinov 1999; Bland- 
ford 2002) and on the numerical procedure described by 
Spitkovsky (2006). The two investigations differ in the outer 
boundary conditions of the computational grid. CKIO imple- 
mented a Perfectly Matched Layer (PML) that guarantees 
absorption/non-reflection of the waves reaching the edge of 
the computational domain. This allowed them to run their 
original 3D pulsar magnetosphere simulation in a small com- 
putational box extending out to only about 2.5 light cylin- 
der radii from the 'star' for as long as necessary to achieve 
steady state. On the other hand, Spitkovsky (2006) and Bai 
& Spitkovsky (2010) did not implement any special outer 
boundary condition but simply used a much larger com- 
putational domain; as such, they followed the evolution of 
the pulsar magnetosphere in the inner 2 light cylinder radii, 
for only about 2 stellar rotations (that is until the waves 
reflected at the boundary reach the distance of 2 light cylin- 
der radii) . Both investigations reached a very similar steady 
state magnetospheric configuration. 

The issues raised in the Introduction suggest the need 
for computing the structure of the magnetosphere on a more 
extended domain than employed so far. Obviously, the im- 
plementation of a PML makes our code more suitable to 
address this problem with the same amount of computa- 
tional resources than any other FFE numerical code cur- 
rently available (see CKIO for details about the code). We 
have run simulations extending out to 10 light cylinder radii 
in all directions for various pulsar inclination angles. Each 
simulation required about 35 hours on a 1000-core super- 
computer. We start off with an inclined magnetostatic dipole 
field at the center of our cartesian computational cube. At 
time t = 0, we set the central dipole in steady rotation. We 
run the simulation for about 2 stellar rotations, which is the 
time needed for the initial transient wave to reach the outer 
boundary of the numerical grid plus another 2 rotations to 
safely reach steady state across our computational grid. The 
resolution of the simulation is 50 grid points per light cylin- 
der radius, i.e. two times higher than that in CKIO (for a 
total of 1000^ grid cells - not taking into account those of 
the PML). We would like to note here that, in distinction to 
our original smaller scale simulations that reached only out 
to 2.5 rzc, our present extended ones include artificial diffu- 
sion. This was required in order to better approach the final 
steady state (we will return to this important point in § 2.3 
and in the end of § 3). More specifically we added artificial 
numerical diffusion terms in the time-dependent Maxwell 
equations 

= -VxE + ryneV'B (1) 

c at 

= VxB J + ryrzcV^E (2) 

c at c 

where 77 is a dimensionless diffusion parameter acting mostly 
in the outer magnetosphere beyond 1.5 ric. The value of 77 
is weighted with |V^B| and |V^E|, and its minimum and 
maximum values throughout the numerical grid are 2 x 10~^ 
and 6 x 10 ~^ respectively. 
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Figure 1. Meridional cross section of the magnetosphere on the poloidal (^i, 17) plane in steady state for various pulsar inclination angles. 
Left panels: meridional electric current density multiplied by (color plot) and direction (arrows). Right panels: electric charge density 
multiplied by (color plot) and meridional magnetic field direction (arrows). Inclination angles from top to bottom: 15°, 30°, 60°, 90°. 
We have considered the mean values of all the plotted quantities over one period in order to reduce the noise level. Light cylinder shown 
with black vertical dashed lines. The black solid lines and the white dashed lines in the right hand column indicate the zero charge lines 
and the corresponding section of the original 'Bogovalov' type current sheet, respectively. 
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Figure 2. Detail of the electric charge density in the poloidal (^i, 1^) plane in steady state for an inclination angle of 30°. Left panel: 
low resolution (25 grid points per light cylinder radius) non-diffusive simulation. Center panel: high resolution (50 grid points per light 
cylinder radius) non-diffusive simulation. Right panel: high resolution (50 grid points per light cylinder radius) simulation with artificial 
diffusion as in Eqs. (1) &: (2). Zero charge (null) surfaces shown with solid black lines. Crossing of null surfaces shown with white line 
elements and circle. Light cylinder shown with dashed lines. 

2.2 Results 

Our results are summarized in Fig. 1 where we plot vari- 
ous physical quantities on the poloidal (/i, ft) plane (i.e. the 
plane that contains both the rotation and magnetic axes) 
in steady state for various pulsar inclination angles. These 
plots were obtained after about 3-4 stellar rotations, beyond 
which an almost steady pattern of currents and fields that 
corotates with the central star is established. Note that in 
Fig. 1 we have considered the mean values of all the plot- 
ted quantities over one period in order to reduce signifi- 
cantly the noise level. On the left, we plot the distribution 
of meridional electric current density multiplied by (color 
plot) and direction (arrows), and on the right the electric 
charge density multiplied also by (color plot) and merid- 
ional magnetic field direction (arrows). The black lines in the 
right hand column denote the zero charge lines. The incli- 
nation angles are 15°, 30°, 60°, and 90° from top to bottom 
respectively. Note that, in FFE, when a corotating steady- 
state is reached, the poloidal electric current and magnetic 
field components should be parallel everywhere. We do ac- 
knowledge that this condition is satisfied to a high degree in 
our simulations but obviously not 100%. The situation was 
much worse in our lower resolution -25 grid points per light 
cylinder radius- simulations. We offer some idea about the 
cause of this effect in the end of § 3. 

The most noticeable element of our solution is that 
an undulating equatorial current sheet forms as predicted 
by Bogovalov (1999, hereafter B99). This is the 3D gen- 
eralization of the equatorial current sheet discovered in 
the axisymmetric investigation of Contopoulos, Kazanas & 
Fendt (1999). The white dashed line in the right hand col- 
umn of Fig. 1 indicate the section of the original 'Bogo- 
valov' type current sheet on the (/i, ft) plane. The current 
sheet survives as far as the outer boundaries of our numerical 
simulation for several stellar rotations and is not destroyed 



by instabilities or reconnection (see however also § 2.3, 3). 
It originates at the positions where the closed line region 
'touches' the light cylinder. The current density at those 
positions is clearly non-uniform. The opening half-angle of 
the current sheet above and below the rotational equator is 
about equal to the pulsar inclination angle, as predicted by 
B99. 



2.3 Numerical or physical instabilities? 

CKIO and iBai Suitkovskvl (|2010l ) have already identified 
regions of enhanced electric charge density double layers 
(positive and negative) described as 'knots' just outside the 
light cylinder. Obviously, the electric field is greatly en- 
hanced between such layers. The results presented in CKIO 
did not include any artificial diffusivity. In the present study 
we realized that this behavior is more prominent in higher 
resolution non-diffusive simulations. In Fig. 2 we plot de- 
tails of the electric charge density distribution in the poloidal 
plane (/i, ft) corresponding to a = 30° as obtained in the low 
(Fig. 2a) and high (Fig. 2b) resolution non-diffusive simu- 
lations, and in the high resolution simulation with artificial 
numerical diffusivity (Fig. 2c). In Fig. 2a we observe that 
the region of strong positive and negative electric charges 
is smoother than in Fig. 2b. As a general trend, in higher 
resolutions numerical instabilities tend to be suppressed, 
whereas physical instabilities tend to manifest themselves. 
The fact that the feature observed in Fig. 2a is enhanced in 
simulations of increased resolution signifies an instability of 
physical rather than numerical underpinnings. Nevertheless, 
we cannot yet rule out that it may also be due to some kind 
of numerical instability. In Fig. 2c we see that the introduc- 
tion of the diffusion term in eqs. (1), (2) clearly smooths out 
this feature. 

We realized that such 'knots' appear where two dif- 
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ferent zero charge surfaces cross each other on the cur- 
rent sheet. In Fig. 2 the sign of the charge density on the 
poloidal plane reveals this strange topology: the crossing of 
two zero charge (null) surfaces creates an X-point on the 
current sheet. Above and below this point along the current 
sheet double layers of opposite electric charges are formed. 
We have also found that dissipative solutions of the mag- 
netosphere above and below this point exhibit electric field 
components parallel and antiparallel to the magnetic field, 
respectively. 

The question with this feature is whether it represents a 
real structure or an artifact of our numerical scheme. Is this 
feature introduced by the non-diffusive scheme (Fig. 2a, b) 
or is it that the diffusion introduced in producing Fig. 2c 
smooths out a really singular feature of the solution? The 
fact that two independent numerical schemes (CKIO; Bai & 
Spitkovsky 2010) exhibit similar behavior and that this is a 
rather extraordinary region (the confluence of 2 zero charge 
surfaces on a current sheet) supports the view that this may 
indeed be a real feature. 



3 THE PULSAR ANTENNA RADIATION 
PATTERN 

Overall, the magnetic field assumes an open split-monopole- 
like structure that corresponds to the radiation field of an 
antenna with E almost perpendicular to B, and both E and 
B dropping like 1/r on average with distance from the cen- 
tral antenna (in our case, the rotating central dipole). We 
remind the reader that in any radiation pattern, a character- 
istic distance exists that separates the inner from the outer 
radiation field region, the so-called Fresnel distance defined 
as 

rpresnel = / X ^ Tic ■ (3) 

Here, ric = cP/ (27r) is the light cylinder radius; / ^ 2ric is 
the linear extent of the 'pulsar antenna'; and A = 27vric is the 
wavelength of the radiation as seen in the electromagnetic 
field patterns of Fig. 1. Our simulation box extends beyond 
a few times the above distance, and therefore, it is suitable 
to study the pulsar antenna far- field radiation pattern. 

While the Poynting flux of the rotating magnetosphere 
averaged over the entire 47r solid angle decreases as ^ 
as expected from the conservation of electromagnetic radia- 
tion, this flux is not distributed uniformly over the sphere. 
In Fig. 3 we plot the Poynting flux distribution on the Moll- 
weide (elliptical) projection of a sphere of radius r = 8ric 
centered on the central star (^ = along the axis of rota- 
tion). One sees that the Poynting flux displays local max- 
ima. Obviously, in steady-state, the radiation pattern should 
corotate with the star. For a = 90° the 50% of the Poynt- 
ing flux (black dotted curve of Fig. 3) is emitted from only 
O.Gtt sr. 

Ardavan (1998) proposed that singular electromagnetic 
radiation features take place in vacuum when a pattern of 
electric charges and currents corotates superluminally (only 
the pattern corotates; nothing is moving superluminally). In 
the Introduction, we acknowledged that this was one of our 
original motivations for extending the simulations of CKIO 
to much larger distances. Preliminary results hinted that the 
local Poynting flux in regions of enhanced electromagnetic 



field discussed in § 2.3 drops as 1/r (and not as 1/r^) with 
distance r, as anticipated in vacuum by Ardavan (1998). 
This was very encouraging! However, the fact that certain 
features of the simulation were noisy and not consistent with 
steady-state led us to increase the spatial resolution by a fac- 
tor of two and to introduce diffusivity in order to stabilize 
the numerical scheme. As a consequence any singular elec- 
tromagnetic radiation feature was killed. 

Unfortunately, we are still not in a position to com- 
pletely dismiss the possibility that the physical solution in- 
cludes regions of enhanced electromagnetic field (radiation 
caustics?) and that the introduction of (albeit small) diffu- 
sion does not artificially suppress their development. Fur- 
thermore, diffusivity introduces also a gradual decrease of 
E with distance below its theoretical steady-state corota- 
tion value {r sin 0/ric) Bp, and therefore, magnetic field lines 
gradually loose corotation with the central star. This man- 
ifests itself in Fig. 1 through the inability of our numerical 
procedure to satisfy 100% the steady-state requirement that 
the poloidal components of the electric current and the mag- 
netic field must be parallel everywhere. 



4 CONCLUSIONS 

In this paper, we presented the results of extensive numer- 
ical simulations of the 3D force-free ideal pulsar magneto- 
sphere out to a distance of about 10 times the light cylinder 
radius. Our results clearly show that the undulating equato- 
rial current sheet survives over that distance and over sev- 
eral stellar rotations, implying that it represents an almost 
steady feature of the magnetosphere. The introduction of 
a small artificial diffusivity was necessary for the stabiliza- 
tion of the numerical scheme. Using this scheme the current 
sheet appears smooth and very similar to the one given an- 
alytically in B99. Overall, it is primarily charged positively 
(negatively) for an aligned (counter- aligned) rotatofl 

Interestingly enough, the regions of strong charge den- 
sity first seen in CKIO are enhanced in higher resolution non- 
diffusive simulations. These regions are formed near the lo- 
cus of the crossing of two zero-charge surfaces. A destabiliz- 
ing behavior more prominent in high resolution simulations 
is observed around these regions. Nevertheless, this behav- 
ior is suppressed after the introduction of artificial diffusiv- 
ity (necessary for the stabilization of the numerical scheme), 
thus raising questions about its true nature (physical or nu- 
merical) . 

We present also the pulsar antenna spiral (for oblique 
rotators) radiation pattern. The Poynting flux is not uni- 
formly distributed over the entire sphere. Its maximum is 
along a spiral direction near the (rotational) equator. The 
Poynting flux corresponding to the smooth high resolution 
solutions produced by the diffusive simulations drops every- 
where as implying that enhanced electromagnetic fields 
similar to those proposed by Ardavan (1998) are not present. 
Their existence, however, cannot yet be conclusively ruled 
out since such effects are expected to be sensitive to diffu- 
sivity. Further investigation with much higher resolution is 
required. 

^ Obviously, the concept of alignment gradually loses its meaning 
as we approach 90° inclination angles. 
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Figure 3. Distribution of the Poynting flux in MoUweide (elliptical) projection of the sphere r = 8ric centered on the central star for 
the inclination angles indicated in the panels. The linear color scale ranges from up to 4 times the corresponding mean Poynting flux 
value. The highest 50% of the values are contained inside the dashed black lines, and the highest 10% inside the solid black lines. The 
maximum of the Poynting flux forms a spiral pattern near the (rotational) equator. Note that the 50% of the Poynting flux in a = 90° 
case (dotted curve) is emitted only through 15% of the entire 47r solid angle. The same area corresponding to the a = 15° case is 40% 
larger. The white dashed lines trace the locus of the current sheet. 
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